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Abstract. We show how matter can be included in a constrained ADM-like 
formulation of the Einstein equations on constant mean curvature surfaces. 
Previous results on the regularity of the equations at future null infinity are 
unaffected by the addition of matter with tracefree energy-momentum tensor. 
Two examples are studied in detail, a conformally coupled scalar field and a Yang- 
Mills field. We first derive the equations under no symmetry assumptions and then 
reduce them to spherical symmetry. Both sectors (gravitational and sphaleron) 
of the spherically symmetric Yang-Mills field are included. We implement this 
scheme numerically in order to study late-time tails of scalar and Yang-Mills 
fields coupled to the Einstein equations. We are able to evolve spacetimes that 
disperse to flat space, accrete onto a given black hole or collapse to a black hole 
from regular initial data. The sphaleron sector of Yang-Mills is found to exhibit 
some nontrivial gauge dynamics. 



1. Introduction 

Most current numerical relativity codes are based on the Cauehy formulation of 
general relativity and evolve spacetime on spacelike slices approaching spacelike 
infinity, truncated at some finite distance. At the resulting artificial timelike 
boundary, boundary conditions must be imposed that yield a well-posed initial- 
boundary value problem and are compatible with the constraint equations. In 
addition, for evolutions of isolated systems, gravitational radiation should pass through 
the boundary without causing spurious refiections, i.e. the boundary conditions should 
be absorbing. Considerable progress has been made recently with the construction 
and implementation of boundary conditions for the Einstein equations (see |Tj for a 
recent review article). A fundamental problem remains, however. In general relativity 
there is no well-defined fiux of gravitational radiation at a finite distance upon which 
absorbing boundary conditions could be based. At best one may appeal to linearised 
theory. Gravitational radiation is only well defined at future null infinity J^^ . Thus 
a far more elegant solution to the outer boundary problem is to include in the 
numerical domain. 

We follow Penrose's approach [5] and apply a conformal transformation to the 
spacetime metric, combined with a eompaetifying coordinate transformation that 
maps an asymptotically flat spacetime to a finite domain. In [3] we developed 
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an ADM-like [J formulation of the vacuum Einstein equations on constant mean 
curvature (CMC) surfaces. These are spacehke but approach future nuU infinity 
instead of spacehke infinity; thus one solves a hyperboloidal initial value problem. 
(Note that hyperboloidal surfaces are not Cauchy surfaces; we only obtain the part 
of spacetime to the future of the initial hyperboloidal surface.) Since the Ricci tensor 
is not conformally invariant, the Einstein equations contain inverse powers of the 
conformal factor that are singular at J^~^ . However, in [3] we showed how the 
formally singular terms in the ADM evolution equations can be evaluated at J^'^ in a 
regular way provided the constraints hold and J^~^ is shear free. In [5] this scheme was 
implemented numerically for axisymmetric spacetimes. Long-term stable evolutions of 
a gravitationally perturbed Schwarzschild black hole were obtained and the Bondi news 
function describing the gravitational radiation emitted by the system was evaluated 
at . 

Before continuing we briefly review other hyperboloidal evolution schemes 
and associated numerical studies. The oldest and, arguably, mathematically best 
understood formulation are the regular conformal field equations due to Friedrich 
[6]. These form a symmetric hyperbolic system of partial differential equations that 
contain the Einstein equations as well as evolution equations for the Weyl curvature 
arising from the Bianchi identities. The equations have the remarkable property that 
they are manifestly regular up to . There have been various attempts at numerical 
evolutions based on these equations (see the review articles [7l[8j[9]). With a view to 
the applications considered in the present paper, we mention in particular the studies 
of spherically symmetric scalar field collapse by Hiibner [TOl [H] . Recently a number 
of formulations have been suggested that are based more directly on the Einstein 
equations, as in our approach. Zenginoglu |12j developed a formulation based on 
generalised harmonic gauge combined with a suitable choice of gauge source functions 
at . Bardeen, Sarbach and Buchman [13] derived a tetrad formulation of the 
Einstein equations on CMC slices. First numerical results on initial data for single 
and binary black holes were presented in [Mj [15] . 

In this paper we return to the constrained ADM formulation on CMC slices 
developed in [5] and extend it to include matter sources. The motivation for this 
derives partly from the fact that we wanted to study the late-time behaviour of 
perturbed black hole spacetimes in the context of the full (rather than linearised) 
Einstein equations, including future null infinity. In [5] we correctly reproduced the 
quasi-normal mode radiation emitted by a perturbed black hole but due to limited 
numerical resolution we were unable to resolve the power-law tail expected at later 
times. Therefore, in order to see if our method is suitable to study these phenomena, 
we decided to take one step back and consider spherically symmetric spacetimes, which 
are computationally less expensive to evolve. Because of Birkhoff 's theorem, matter 
is needed in order to have nontrivial dynamics in spherical symmetry. How to include 
matter in a hyperboloidal Einstein evolution scheme is an interesting problem in its 
own right. 

At late times matter fields as well as gravitational perturbations on flat space 
and black hole spacetimes typically decay polynomially in time, a phenomenon often 
referred to as Price's law [l^. This power-law tail is caused by the backscatter off 
the curved background spacetime and/or by the nonlinearity of the matter fields. It 
plays an important role in trying to prove stability of black hole spacetimes and the 
cosmic censorship conjecture [17]. At J^+the fields generally decay at a slower rate 
than at any finite distance, although the closer an observer is to J^'^ , the longer the 
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measured decay rate stays close to the value corresponding to before it ultimately 
approaches the faster finite-distance decay. It can be argued [H] that the decay rate 
at is the relevant one for observers in the astrophysical zone [TH] , in which the 
distance to the source of radiation is very large compared to the time during which 
the signal is observable. 

For scalar fields the tail decay rates can be predicted from linear perturbation 
theory [20( I21j and have been confirmed numerically many times. We mention two 
recent studies that both include in the numerical domain. Piirrer, Husa and 
Aichelburg [TB| evolved the spherically symmetric Einstein-scalar field system in Bondi 
coordinates and determined the power-law tails in subcritical evolutions that disperse 
to flat space. Bondi coordinates cannot penetrate horizons so the decay of the field at 
the horizon of a black hole could not be studied. Zenginoglu evolved the spherically 
symmetric scalar wave equation on a fixed background spacetime (taken to be either 
Minkowski or Schwarzschild), the test field approximation. He used a hyperboloidal 
foliation of spacetime that covers part of the black hole interior as well. In the present 
paper we evolve the coupled Einstein-scalar field equations on hyperboloidal slices 
reaching out to J^~^ . In particular, we are able to study the decay of the field at the 
horizon of a black hole formed in gravitational collapse. 

For Yang-Mills fields the prediction from linear perturbation theory turns out to 
be incorrect: the nonlinearity of the field causes a slower decay [52]. Similarly to the 
scalar field, the Einstein- Yang-Mills system was evolved in Bondi coordinates in [23] 
and in the test field approximation on hyperboloidal slices in [24]. We shall compare 
our results with those studies. 

This paper is organised as follows. In section [2] we extend our hyperboloidal 
Einstein evolution scheme [3] to include matter sources and we re-examine the question 
of regularity at . Two matter models are studied in detail in section [H] a 
conformally coupled scalar field and a Yang-Mills field. In section |4] we reduce our 
formulation to spherical symmetry. The numerical implementation of this system is 
described in section [5] and our results on power-law tails are presented in section \6\ 
Finally we conclude and discuss some directions for future work in section [7] Some 
useful identities for conformal transformations and 3 -I- I decompositions are collected 
in [Appendix A 

2. General formalism 

In this section we consider general matter subject to the condition that its energy- 
momentum tensor be tracefree. We extend the conformal ADM formulation of the 
Einstein equations derived in [5] to include the corresponding matter source terms. 
The question of regularity at future null infinity is re-examined. 

2.1. Matter in the conformal setting 

We use the notation and conventions of [3] . The spacetime metric gfj^i, is written as 

W.9,,^f^-2W7,., (1) 

where '■^-'7^1/ is the conformal spacetime metric and fl the conformal factor. Consider 
now matter given by an energy- momentum tensor Tf^^. We introduce a conformally 
rescaled energy-momentum tensor T^^, via 
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The energy-momentum conservation equations transform as |25j 

where ^''^V denotes the covariant derivative of '^^^g and ^''^V the covariant derivative 
of (^^7. The first term on the right-hand side of ([3]) vanishes by energy-momentum 
conservation. If the energy-momentum tensor is tracefree, 

= 0, (4) 

then the second term on the right-hand side of (jS]) also vanishes and the energy- 
momentum conservation equations reduce to 

(4)^A.. (4)v^f,p = 0. (5) 

Assuming that T^i, is itself regular at J^"*" , which it is if its basic fields satisfy 
conformally regular field equations (as in the matter models considered in section 
O, then (O is manifestly regular at future null infinity . For this reason we will 
restrict ourselves to matter with tracefree energy-momentum tensor. 

Examples of matter models satisfying this condition include the conformally 
coupled scalar field (section lXTj) . Maxwell and more generally Yang-Mills fields (section 
13. 2p and the radiation fiuid, i.e. a perfect fluid with equation of state p = -^p (see also 
[26]). There are of course many matter models that do not satisfy ([4]). Typically 
however, these are considered to be less radiative. For example, a massive scalar field 
is known to fall off faster than any power of radius towards J^"*" [27]. If the matter 
remains bounded away from .y^ , as is expected in many situations of astrophysical 
interest, then of course there is no harm in using the singular matter field equations 
away from J^^ . 



2.2. Conformal ADM reduction of the Einstein equations with matter sources 

As in we decompose the physical and conformal spacetime metric in 3 + 1 form, 

(4)g ^ ~N^dt^ + gij{dx' + X'dt){dx^ + X^dt), (6) 

('')7 ^ -N^dt^ + -fij{dx' + X'dt){dx^ + X^dt), (7) 

where the physical and conformal lapse are related via N = n~'^N, the physical and 
conformal spatial metric via g = ri~^7, and X^ serves both as physical and conformal 
shift. The unit timelike normals to the t = const hypersurfaces in physical and 
conformal spacetime are related by — fln'^. 
The extrinsic curvature is defined as 

Kij = -\Cngij. (8) 

Our condition on the time slices is that their mean extrinsic curvature be constant, 

g'-^K.j = -K = const (9) 

with K > so that the slices approach future null infinity (hence our slightly awkward 
sign convention in ©). Wc choose to work with the traceless part of the ADM 
momentum, which is given in terms of the extrinsic curvature by 

TT*"^^ =-M,(gV'-kV)A'^i, (10) 



where p,g = det{gij) . For later use we also define the conformal extrinsic curvature 

Cij ~ -^Cnjij. (11) 
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Its trace C = ^^■'Cij is related to the variable F — -~2NC introduced in [3]. 

Following [28) we define the following projections of the energy-momentum tensor, 

p = n'^n%,, J' = -g'^^n''T^,, S,, = g.^g^-'T^,, (12) 

and similarly for the conformally rescaled energy- momentum tensor ([2|), 

p = fL'^h''f^,, r = ^f'^h''f^,, = 7,''7/f^,. (13) 

Clearly these quantities are related by 

p = n^p, r = n'^r, s^J = n^s.j. (m) 

The Einstein equations are 

Gi_,„ = KTfj^^ (15) 

with K = Stt in geometric units (i.e. Newton's constant and the speed of light 
G = c = 1). They split into evolution equations and constraints. The evolutions 
equations are (cf. equations (2.10), (2.21) and (2.22) in [3]) 

Cnn = ~^{K + nC), (16) 
£r^l^J = 2p~^-f,kjji7r''''' - |7,,C, (17) 



-\-p. 



(18) 



Here and in the following, indices on quantities carrying a tilde are to be raised and 
lowered with 7, tr denotes the traccless part w.r.t. 7, Rij is the Ricci tensor of 7, and 
C denotes the Lie derivative. Note that tt''^*-' is a tensor density of weight one and 
hence 

Cnir'"^ = N-^ [dtn'"^ - (X'^7r*"^Xfc + X\k7r''''^ + X^fcTT*"'^'] . (19) 
The Hamiltonian and momentum constraints arc (cf. equations (2.9) and (2.7) in [3]) 

= -inv^Vjfi + 67'Jrj,,i7j - n^R - ^k^ + n^n-^j.,kijiTT^"^Tr^'''^ + 2kVl^p, (20) 

Q = Vj{n-^-K^"'^) + np^J\ (21) 

Preservation of the CMC condition, dtK = 0, implies an elliptic equation for the 
conformal lapse (cf. equation (2.13) in [3]), 

+ lNfl^fi-^j,kljin''''iTr'''^' + ^kM^S + 2p) = 0, (22) 
where S = "f^^ Sij. 

We also need to specify the spatial coordinates. In [3] we imposed a spatially 
harmonic gauge condition, which yielded an elliptic system for the shift (equation 
(2.15) in [3]). However, other choices are possible. For example, in [5] and in section 
m below we use coordinates adapted to a spacetime symmetry. 

There is a residual conformal gauge freedom inherent in the decomposition ([T]). 
In [3] we fixed this by requiring the conformal scalar curvature R to be constant; this 
resulted in an elliptic equation for F = —2NG (equation (2.12) in [3]). In section |4] 
the spherically symmetric conformal metric will be taken to be flat; this eliminates 
the conformal gauge freedom. 
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2.3. Regularity at future null infinity 

The evolution equation for the traceless momentum (jl8p and the elhptic equations 
(|20p - (P^ contain inverse powers of the conformal factor that are singular at ^+ . 
However in [3] we showed, for the vacuum case, that the formally singular terms in 
(|18p can in fact be evaluated at .y^ in terms of regular conformal quantities, provided 
the constraints are satisfied and ^+ is shear free (see also (25). We also showed that 
these regularity conditions are preserved under the time evolution. 

Our analysis exploited the fact that the elliptic equations are degenerate 

at ^+ . This allows one to evaluate the first few radial derivatives of the fields 
explicitly at J^"*" , where radius r refers to a coordinate on a given spatial slice that is 
constant on the cut of the slice with . In particular, we obtained expressions for 
the first three radial derivatives of J7, the first two radial derivatives of N and the first 
radial derivatives of tt*'' at J^^ . This information was sufficient in order to evaluate 
the formally singular terms in (jlSp by applying L'Hospital's rule. 

It is easy to see that the addition of matter does not affect those results, essentially 
because the matter terms in ((20 |) -([22 |) are multiplied by sufficiently high powers of 
that vanish at ^+ even after a certain number of derivatives are taken. More 
specifically, the third radial derivative of ft at J^^ was obtained in [3] by evaluating 
the Laplacian of (|20p . Now the Laplacian of the matter contribution to ([201 is 0{ft^) 
and hence does not contribute at =y+ . Clearly, the matter contributions to lower 
derivatives of are multiplied by even higher powers of and vanish at as well. 
The first radial derivative of tt^'^" was obtained by first multiplying (|2T|) by ft^ and 
then taking a radial derivative. After this operation the matter contribution is 0{fP) 
and hence, again, does not contribute at . We found the second radial derivative 
of N by forming the linear combination x ([22)) + jfl~^N x ([20]) and then taking 
a radial derivative. The matter contribution to the resulting expression is 0(17^) and 
hence does not contribute at . Thus all the expressions we derived in order to 
evaluate the formally singular terms in ([181 at are unchanged by the addition of 
matter. The matter term in ([T5)) itself is 0(r2^) and thus vanishes at . 

3. Examples of matter models 

In this section we consider two matter models in detail, a conformally coupled scalar 
field and Yang-Mills theory. We write the matter equation of motion in 3 + 1 form 
and work out the projections of the energy-momentum tensor that appear as source 
terms in the 3 + 1 form of the Einstein equations. 

3.1. The conformally coupled scalar field 

The action for the Einstein-scalar field model with conformal coupling is given by 



This differs from minimal coupling by the additional last term containing the four- 
dimensional Ricci scalar. As we shall see shortly, the advantage of the conformally 
coupled model is that it yields a conformally invariant evolution equation for the scalar 
field. We note however that the conformally coupled Einstein-scalar field equations 
are equivalent to the minimally coupled ones in the following sense |30[ llOj : Suppose 




(23) 
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{(ff^, (4)(^cc^j ^ solution to the conformally coupled equations. Then {(j)™'^^ '"^•'s^J^) is 
a solution to the minimally coupled equations, where 

= V6arctanh(ir^), W^„.^c ^ _ i(^cc)2] (4)^cc ^ (34) 

Varying (j23p w.r.t. (p yields the equation of motion 

□0-i(4)i?0 = O. (25) 

Defining a conformally rescaled field 4> = ^^id using the transformation of the 

Ricci scalar (|A.4p . we find 

i = r2-3(n0_ 1 (4)i^^) ^ 0. (26) 

The left-hand side of this equation is manifestly regular at J'^ . 

Varying (j23p w.r.t. ^"^^ g^u produces Einstein's equations G^u = KT^,y with energy- 
momentum tensor 

T^u = - \4> W ^ 1^2 (4)^^^ _ 1 (4)^^^( (4)gP-</,_^0^, + 102 (4)^)^ (27) 

where an overall factor of | has been absorbed into the definition of 0, and the equation 
of motion has been used to eliminate a 00 term. This can also be written as 
T^y = n^Tfj^i^, where the conformally rescaled energy-momentum tensor has exactly 
the same form, 

Using the 3 + 1 decomposition of the scalar field Hessian ( [Appendix A. 3 1 and the 
conformal spacetime Ricci tensor ( [Appendix A.4 ), the equation of motion (|26p can be 
written as 

= -£l<p + V,V> + CCn4> + N~^f^N^,4>^j 

-i0 (^--2N-^V'V,N + R + ^C^ + fi~^j,kljin'''''n'''^' - 2/:sC)(29) 

(Recall that we are raising and lowering spatial indices with 7jj .) In general, the time 
derivative of the conformal mean curvature C cannot be evaluated explicitly during 
an evolution, e.g. if F = —2NC is determined by solving an elliptic equation as in [5]. 
This problem can be avoided by introducing a new variable tp and writing (j29p as a 
first-order (in time) system 

UP = + \CP, (30) 

(-27V-iV'V,iV + R + fi-^-f^kljiTT'^'TT'^"') . (31) 
The matter source terms in the geometry equations evaluate to 

p = -i0V,V^0+ i0,;0^' + + i02^_ l02^-2^^^,^^.^^tr,j^trfcZ^ (32) 

J' = Y'i-hi' + - (33) 

S^P, (34) 

+ \pp {p-^C-y"^ + 2^-27H7r'"'^7r''-J' - iV-i [V'V^TV]*^ 7?"^^) . (35) 

The equation of motion ([^n|) has been used in order to eliminate a £?0 term in ([5^ . 
Note that (|55|) contains £ft7r''^*-' . When (|55|) is substituted in the evolution equation 
(fT8|) for the traceless momentum, the resulting equation must be solved for LnW^'^'^^ . 
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3.2. Yang- Mills theory 

The fundamental field of Yang-Mills theory is a connection or vector potential 
which in addition to the spacetime index /i carries an index (a) referring to the internal 
gauge group. We use internal indices from the beginning of the Latin alphabet a, 6, . . . 
ranging over 1, 2, . . . , A^, where N is the dimension of the associated Lie algebra, e.g. 
TV = 3 if the gauge group is SU{2). Repeated internal indices arc summed over. 

Since the Yang-Mills equations are conformally invariant, we choose to work 
directly in the conformal spacetime here. The Yang-Mills connection and its conformal 
counterpart are identical, A|f ^ = A^"-* . We regard 

= (36) 

as a gauge variable that can be freely chosen. 

The Yang-Mills field strength tensor is defined as 

F^^) = d^Ai^^ - + /""M^if^) , (37) 

where the symbol /"'"^ is totally antisymmetric in the chosen basis of the Lie algebra. 
For SU{2) we may write /"'"^ = g[abc], where g is a constant; we set g = —2 for 
our numerical evolutions in section [6] The symbol [abc] is totally antisymmetric with 
[123] = I. 

We introduce the electric field 



and magnetic field 

B^^'^^ = ^^,,e^^'^B'^l^='^[^jk]F^^\ (39) 

where 

B^^^F^f (40) 

denotes the spatial part of the field strength tensor and Cijk is the alternating symbol 
associated with the conformal spatial metric -fij. Note that 2?* and arc vector 
densities of weight one. 

The energy-momentum tensor is given by 

T^.. = Fl^;^F,P^-^ - i ^^hi^.H^^FP-^''^. (41) 

The projections defined in ([T3| take the form 

p = \ii-^{V'^''^'D'f^ + ^'(°)^|")), (42) 

J, (43) 

5 = p, (44) 

Energy-momentum conservation yields the Yang-Mills field equations 
Contracting this with n^, we obtain the constraint 

a^P'-l-^) + fabcJ^{b)f,^(c) ^ 0; (47) 
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contracting with 7,^' instead, 

An evolution equation for the vector potential is obtained from the definition of Fq""^ , 

Equations and constitute the evolution equations for the independent 

variables Ji,-"'' and Note that i?,-,"-' is to be expressed in terms of Jl,-"' and 

■ ,1' _ . ^ 

its spatial derivatives using definitions dST]) and 
4. Reduction to spherical symmetry 

In this section we reduce the formulation of the Einstein equations presented in section 
12.21 and the specific matter models of section [3] to spherical symmetry. Special care is 
needed in the case of Yang-Mills theory. 

Einstein equations in isotropic gauge 

We work in spherical polar coordinates t,r,9,tp. All fields are functions of t and r 
only. Partial derivatives w.r.t. time and radius are denoted hy'^dt and ' = dr- 
ill spherical symmetry we may take the spatial conformal metric to be fiat, 

7y =diag(l,r2,r2sin2 0). (50) 

The shift vector has a radial component only, and we set X^' = rX. The traceless 
momentum is diagonal and has only one independent component, 

^trrr ^ _2r^T,t,ee ^ -2r'^smHTr'"^^ = r^sm0TT. (51) 

(We have pulled out a factor of /x^ = sin 6 in the definition of the variable tt in order 
to avoid frequent divisions by /x^ in the final equations.) 

Preservation of ([50]) under the evolution equation pT]) implies the spatial gauge 
condition 

rX' = -f^TT. (52) 
This implies that tt is actually 0{r'^) near the origin and so we set 

TT = r^TT. (53) 
We also note the expression for the conformal mean curvature in our gauge, 

C = N^^ViX' = -^r^n + 3N-^X. (54) 

The reduction of the Einstein equations is as follows. The Hamiltonian constraint 
((20|) reads 

- Ann" + en'^ - 8nr~^n' + foV^Tr^ - ^k^ + 2Kn*p = o. (55) 

The momentum constraint (|2ip is 

n{rTT' + 5tt) -2rn'TT + Kn^r-\r =0. (56) 
The CMC slicing condition ([22]) is 



n'^N" + snn'N' - 2n^r-^N' - ln''^N + Ink"^ + ^■Nn'^r^-n'^ 

+ \KNn^{S + 2p) (57) 
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The evolution equation for the conformal factor (IT51) is 

n = rXn' -Xn + N{^nr^T: - ^K) (58) 
and that for the tracclcss momentum ^TE\i 
TT = rXn' + 3Xtt + ^r~^{r~'^N'y 



-N 



4..1.1. Regularity at the origin. The fields $7, tt, N and X are even functions of r. 
It follows from the general behaviour near the origin of spherically symmetric tensor 
fields [31] [32] that the r-component of a vector J* is 0{r) and the rr-component 
of a tracefree symmetric tensor S'*'''-' is 0{r'^) near the origin; this can also be seen 
explicitly for the specific matter models in the following subsections. Hence the terms 
p, r~^J^, S and r~^S^'^^'^ appearing as sources in ([55|) -(|59 | are also regular even 
functions of r. Keeping this in mind, equations ()55|) -()59 p have been written in a form 
that is manifestly regular at the origin. 

4-1-2. Regularity at future null infinity- For completeness and for later use we also 
state the results of our regularity analysis at here. In the following = denotes 
equality at J^^ . By definition, 

= 0. (60) 
The Hamiltonian constraint (|55p and momentum constraint (j56p yield 

O' = rO" = - \K, n'" = 0, (61) 

TT = tt' = 0. (62) 
From the CMC condition ([57)1 we obtain 

rN' ^ r^N" ^ N. (63) 

Note that the value of N at ^+ can be freely chosen; see section [5^ for our particular 
choice. Preservation of (|60p under the evolution equation (j58p implies 

rX ^ -N. (64) 

The isotropic gauge condition (|52|) further yields 

X' = X" = 0. (65) 

Given these results it can be verified using L 'Hospital's rule that the evolution equation 
(|59p for the traceless momentum reduces to 

dt^ - 0, (66) 

consistent with 
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4.2. The conformally coupled scalar field 

The scalar field evolution equations pop - (PT|) take the form 

4>:^rX4)' + X4> + N{~l4>r^7r + i;), (67) 

■4, = {rX^ + N4>'Y + X^ + N{2r~^4)' - \4>r'^TT^ + i^ir^Tr) + \4>{N" + 2r-^N'). (68) 

The source terms (I32l) - (j35p appearing in the Einstein equations are 

~p = -i0(0" + 2r-i0') + + - (69) 

r-^r = -r-^4>'iP + i0r-V' - ^4>4>'r7T - \4>'^{rTi' + Stt), (70) 

~S = P. (71) 



2 ■ 

7V-i(7T - rXn' - 3X7r) + ir^Tr^ - |iV-V-i(r-i7V')'j . (72) 

Notice again that these equations are manifestly regular in the origin, given that 
(f) and -0 are even functions of r. 

4-3. Yang-Mills theory 

In the case of Yang-Mills it is more convenient to work in Cartesian coordinates, 
so indices ... will refer to Cartesian coordinates in this subsection. Hence the 
conformal metric is now 7^ — 5ij, with fij — 1^ and the shift vector is X^ = Xx^. 

We take the gauge group to be SU (2). The most general ansatz for the spherically 
symmetric Yang-Mills connection is of the form [331 1311 ISS] 

i**") = [aij]x^F + [x^x' - r^S''')H + S^'L, = Gx". (73) 

Here F, H, L and G are functions of t and r only. 
Similarly, we write the electric field as 

V'i'^) = {ai2\x^BF + {x^x' - r'^5''^)DH + S'^'Dl- (74) 

The evolution equations for the vector potential (|49p read 

F = rXF' - NDf + 2XF + g{G - XL){L - r^H), (75) 
H = rXH' - NDh + r-^G' + g{G - XL)F + X{2.H - r-^L'), (76) 
L = - NDl + rG' + G. (77) 

The evolution equations for the electric field (|48)) are 

Df = {rXDF - NF'Y + 2XDf + g{DL - r^DH){G - XL) 

+N[-Ar-^F' - 2gLrH' + gr-^L'{3L - r^H)] 
+r-^N'[-2F + gL{L - r^H)] + Ng[-3F^ - r^H^ - AHL 

+g{r^F^ - 2r^FHL + r^FH^ + 2FL^)], (78) 

Dh = {rXDn - NH')' - r^^XDL - Nr^^L')' + r"^iV'(-3iJ + gFL) 

+gDF{G + XL- 2Xr^H) + XDh{1 + 2gr'^F) - 2gXFDL 

+N[-'ir-^H' + g{2HrF' - 2FrH' + 'iFr-^L')] 

+Ng{-4:FH + g[F'^iL + r^H) + r^H^ - r^H'^L]}, (79) 

Dl = X[2il - gr^F){DL - r^Du) + 2gr^DF{L - r^^H)] 

+2N{1 - gr^F){3H + rH' - r-^L' - gFL) 

-2Ng{L - r^H)[rF' - gr^H{L ~ r^H) + 2F - g{L - r^H)'^]. (80) 
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The Yang-Mills constraint (gT]) is 

r-^D'^ + 2gF{DL - r^Dn) + 2Dh + 2g{r^H - L)Df = 0. (81) 
The components of the magnetic field 

= [aij]x^Bp + {x'^x' - r^S"')BH + S^'Bl (82) 

evaluate to 

Bp = -3H - rH' + r-^L' + gFL, (83) 
Bh =r-^F' - gH{L-r^H)+gF'^, (84) 
Bl ^ -2F + g{L - r^Hf + gr'^F'^. (85) 

In terms of these, the matter sources in the Einstein equations are given by 

p = S^ ^[3Dl - 2r\2DLDH - - r^D%) 

+3Bl-2r^2BLBH - Bl-r^B%)], (86) 

r-^r = 2[DlBf - DfBl + r^iDpBH - DhBf)], (87) 

~-2^tr.r ^ _2^2DlDh - 0% - r^ Djj + 2BlBh -Bl- r^B%). (88) 

We are free to impose one gauge condition on the Yang-Mills connection, and we 
choose radial gauge L = 0. The evolution equation (j77p for L now turns into an ODE 
determining G, 

rC + G- NDl = 0. (89) 

If the initial data are such that H = Dh = then Dl = is a solution 
to ([8T|) . and G ~ (i.e. temporal gauge) is a solution to (f89| . (There are other 
solutions to these last two equations but the boundary conditions we impose will 
single out the given ones.) These conditions are preserved under the time evolution, 
i.e. H = G ^ Dh = Dl = at all times. This reduced system is sometimes 
referred to as purely magnetic or as the gravitational sector of the Einstein- Yang-Mills 
equations, whereas the remaining system is called the sphaleron sector |36j . While 
the gravitational sector has been widely studied numerically (e.g. in [371 [53] ) , we 
are not aware of any numerical simulations using the full system. 

We also note that the ansatz ([75)1 is often written (e.g. in [35]) in a different 
gauge known as Abelian gauge, 

A = WT^ de + {cot e t'' +wt'^) sine d9, (90) 

where denote the generators of SU{2). The two ansatze are related by a (singular) 
SU{2) gauge transformation (see e.g. [38]). 



5. Numerical method 



In this section we describe the numerical methods we use in order to solve the 
spherically symmetric reduction of the Einstein equations on hyperboloidal slices 
derived in section[4j We begin by summarising our evolution scheme. Particular care is 
spent on the boundary conditions. We then describe the finite-difference discretisation, 
time integration method and elliptic solver. Finally we explain how apparent horizons 
are detected and excised. The code has been implemented in Python. 



Hyperboloidal Einstein-matter evolution and tails for scalar and Yang-Mills fields 13 



5.1. Evolution scheme 

Our fundamental geometric variables are f2, tt, N and X. The matter variables are 
either and ip for the scalar field, or F, H, G, Dp, Dh and Dl for Yang-Mills. 
The numerical evolution proceeds as follows. At each time step, we solve 

(i) the Yang- Mills constraint (gT]) ior D^, 

(ii) the Einstein constraints ([55]) and ([55)1 for O and tt. If the source of the 
momentum constraint is independent of tt, as is the case for Yang-Mills, then the 
two equations may be decoupled by setting tt = fi^P, which turns the momentum 
constraint ([55)) into 

rP' + 5P + Kr-\r ^ 0. (91) 

This is solved first and the solution for P is then substituted in the Hamiltonian 
constraint (|55p . For the conformally coupled scalar field, the two constraints 
cannot be decoupled in this way because depends on tt (cf. equation (|70)) V In 
this case the Einstein constraints must be solved as a coupled system. 

(iii) the slicing condition ([57]) for N, 

(iv) the isotropic spatial gauge condition (|52p for X, 

(v) the Yang-Mills radial gauge condition ([89]) for G. 

(Obviously steps (i) and (v) are only included for Yang-Mills matter.) The matter 
variables are then evolved to the next time step using either ([67)) and ([68)) for and 
ijj or ([75]), ([76]), ([78]) and ([79]) for F, H, Dp and Dh- 

The evolution equations ([58]) . ([59| and ()80p for fi, tt and Z?l are not solved 
explicitly but are evaluated during the numerical evolution as a consistency check. 

5.2. Boundary conditions 

We discuss the boundary conditions for each of the equations in turn. We need to 
distinguish between two cases, a regular centre and a black hole with an inner excision 
boundary. The outer boundary is always placed at . 

The Yang-Mills constraint ()8ip is a first-order ODE for Dp that requires one 
Dirichlet boundary condition. The value of Dp at either the origin or the excision 
boundary is obtained by evolving Dp there according to its evolution equation ([80]) . 

The Hamiltonian constraint (|55l) is a second-order ODE for fl, and we regard 
it as a two-point boundary value problem. For a regular centre, the inner boundary 
condition follows from the fact that ft is an even function of r: 51' = at r = 0. For 
an excised black hole, we use the evolution equation ([55| to find the value of ft at the 
excision boundary. The outer boundary condition at is = in both cases. 

The momentum constraint ([56| is a first-order ODE for tt. At an excision 
boundary we use the evolution equation ()59p to obtain a Dirichlet boundary condition 
for TT there. For a regular centre, however, there is a unique solution to (|56p that is 
regular at r = 0. This can be seen immediately by evaluating ([55]) at r = 0, where the 
equation implies a Dirichlet condition for tt. Thus no additional boundary condition 
is imposed at r = 0. 

The CMC slicing condition ([57]) is a second-order ODE for N that we regard as a 
two-point boundary value problem. At a regular origin ?■ the boundary condition 
is N' = 0. At an excision boundary we choose to freeze the value of N from the time 



Hyperboloidal Einstein-matter evolution and tails for scalar and Yang-Mills fields 14 



when the black hole is excised. At ^+ we choose the conformal lapse such that time 
t corresponds to proper time r of an observer at . This is given by 

- dr^ = n-^{-N^ + r^X^)dt^. (92) 

Using the results of our regularity analysis (section 14. 1.2[) . the formally singular term 
on the right-hand side can be shown to have a regular limit at J^"*" , 

n-^{N^ ~r^X^) =9r-^K-^N^. (93) 

Hence choosing t to coincide with r at ^+ corresponds to setting 

N = ^Kr. (94) 

The isotropic spatial gauge condition ([5^ is a first-order ODE for X. The 
boundary condition for X at ^+ follows from preservation of = under the 
evolution equation (|58| . namely, X ^ — r^^N. 

Finally we consider the Yang-Mills radial gauge condition (|89|) . This first-order 
ODE for G has a unique solution that is regular at r = 0. (Evaluating the equation 
at r = implies a Dirichlet condition for G.) Thus no additional boundary condition 
must be imposed at r = 0. For an inner excision boundary we choose to freeze the 
value of G there from the time of excision. 

Since the outer boundary at is null and the inner excision boundary is 
spacelike, the evolution equations (in particular those for the matter variables) do not 
require any boundary conditions there. Near a regular origin the evolution equations 
need not be modified either. The way we discretise the equations near the boundaries 
is described in the following subsection. 



5.3. Discretisation 

The numerical domain is an interval r G [0, 1] for a regular centre or [r„ii,i, 1] for an 
excised centre. In both cases the outer boundary at r = 1 corresponds to J^^ . 

For solutions containing a black hole, the fields typically have steep gradients close 
to the inner (excision) boundary and hence it is advisable to introduce a non-uniform 
grid. We do this by introducing a new radial coordinate x with respect to which the 
numerical grid is uniform, combined with a map 

r : [0,1] ^ Ki„,l], a;^r(x), (95) 

that is steeper near x ~ 1 than near a: = 0, thus providing more resolution close to 
the inner boundary. As in [5] we choose 

r{x) = Qix"^ + (1 - '^min - Qi)x + r,ni„, (96) 

where ^ Qi < 1 is a constant, typically Qi = 0.5. For evolutions with a regular 
centre, we use r(x) ~ x for the numerical results shown here. More generally, we have 
implemented a map 

r{x) ^ Q2X + {1 ~ Q2)x\ (97) 

with ^ Q2 1 a constant. Note that since r{x) in (|97)) is an odd function of x, 
the parity of a grid function with respect to x will be the same as with respect to r. 
(This is not the case for (|96|) but parity considerations are only relevant for regularity 
at the origin.) 

The interval x € [0, 1] is covered by equidistant grid points. For a regular centre 
we use a grid that is staggered about the origin, 

x^ = ii+^)h, O^is^N, (98) 
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whereas for an excised centre we use an unstaggered grid 

X, = ih, Oi^i^N, (99) 

so that ro = r(xQ) — r{0) ~ r^in hes on the excision boundary. In both cases, the 
grid is chosen such that rjv = r{xN) = 1- 

Derivatives with respect to x are discretised using fourth-order accurate finite 
differences. Exphcit expressions for the finite-difference operators can be found in 
Appendix C of [S] . We use one-sided differences near ^+ and the inner excision 
boundary. For a regular centre, the usual centred finite-difference operators are used 
near the origin by formally extending the staggered grid (|M|) to negative values of i 
and using the fact that all the evolved variables are even functions of r (and hence of 
x), i.e. replacing 

X-i ^ Xo, X-2 ^ Xi, X-3=X2, ... . (100) 

The reader will have noticed that we have written the highest derivatives in 
the matter evolution equations and ([75)) -([7 ^ in flux-conservative form. When 
discretising these equations, we apply the finite-difference operators in prescisely the 
order in which the terms in the continuum equations are written. This has been found 
to be essential for the stability of the method. 

5.4- Time integration 

A fourth-order Runge-Kutta method is used in order to integrate the evolution 
equations forward in time. At each full timestep, the elliptic equations are solved 
as described below in section 15.51 At the substeps of the Runge-Kutta algorithm, 
we have found it sufficient to extrapolate the solution to the elliptic equations from 
previous timesteps (using a cubic polynomial approximation) instead of solving the 
equations explicitly. Kreiss-Oliger dissipation j39j is added in order to ensure stability; 
see Appendix C of [5]. For a regular centre, the dissipation operator is evaluated up to 
the innermost grid point using the rule (|100p . No dissipation is added at the outermost 
two grid points near and the excision boundary. 

5.5. Elliptic solver 

The elliptic equations listed in section 15. 1[ which are in fact ODEs in spherical 
symmetry, are solved using a combination of a Newton- Raphson iteration and a direct 
band-diagonal solver. The Newton-Raphson iteration is actually only needed for the 
Hamiltonian constraint (|55|) . all the other equations are linear. The matrices arising 
from our fourth-order finite-difference discretisation are pentadiagonal, except near 
and the excision boundary, where one-sided differences with wider stencils are 
used; a few Gaussian eliminations are applied by hand in order to reduce the matrix to 
pentadiagonal form there. We use the Python routine numpy . linalg. solve Jsanded 
(a straightforward generalisation of the Thomas algorithm) in order to solve this 
pentadiagonal linear system. 

5. 6. Apparent horizon finder, black hole excision and Bondi mass 

We detect the formation of an apparent horizon by tracking the optical scalars 6± 
during the evolution. These are defined as follows. Let = Q,^^ N {dt) ^ denote the 
(physical) unit timelike normal to a given t = const slice and — il{d/dr)'^ the unit 
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outward normal to an r = const surface within that shcc. The outgoing and ingoing 
null normals to this surface are 1^. = n^^ zL . The optical scalars or null expansions 
are now given by 

9± = (Ini?).^/^ = iA' - ^nr^Tr±{r-^n - n'), (101) 

where R = r/Q denotes (physical) areal radius. 

An apparent horizon is given by the outermost radius r at which = 0. A 
zero of this function is detected using a standard root finding algorithm (the Python 
routine scipy. optimize. brentq). Even before a zero forms, we can detect and track 
a minimum of 6^ in order to have a better initial guess at where the zero is about to 
occur. 

Once an apparent horizon has formed at a radius r = 7'ah (with corresponding 
mass A/ah — ^Rah), we excise just inside it, at r = 0.9 tah = ?'min- The numerical 
solution is interpolated (using cubic interpolation) to a new grid with inner boundary 
at r — rniin, discarding the part of the solution at smaller values of r. The evolution 
is then continued on the new grid. 

The optical scalars ()101|) are closely related to the Hawking mass 

Mii = ^R{1 + R^O+e^). (102) 

Its limit at is the Bondi mass Mb , which will be a useful quantity to evaluate 
during our numerical evolutions. Although formally singular, the results of our 
regularity analysis (section I4.1.2p imply that the limit at of (jl02p is 

Mb = MuU+ = -^K-^{Kr'TT" + 2r^n'^^y). (103) 

Evaluating the fourth derivative of the conformal factor at J^"*" is prone to numerical 
error; we have found it to be more accurate to extrapolate Mh from the interior. 



6. Numerical results 



6.1. Initial data 

We start the evolution from initial data that are close to either Minkowski or 
Schwarzschild spacetime. The geometry variables {ft, tt, N and X) are first set 
according to the respective vacuum solution, then initial data for the matter fields 
are specified, and finally the constraints and elliptic gauge conditions are re-solved for 
the geometry variables. 

Schwarzschild spacetime in CMC coordinates is given by [101 SI] 

^2 _ ^-^ 2Af ^ 1 2a 

where 



ds^ = -[l- dt^ + —df^-Ydtdf+f^{d9^+sm^0d<j)^),{lO4) 



2M ^V^"^ , , Kf C 



fir)^[l- — + a'j , a(f) = — --, (105) 

and M (mass), K (mean curvature) and C are constants. The radial coordinate f 
needs to be transformed to a new radial coordinate r such that the spatial metric is 
conformally fiat in the new coordinates. This yields the ODE 

dr r 
dr rjir) 
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Since the physical radial coordinate f has infinite range, it is more convenient to work 
with s = 1/f and set 

A{s) = \K - Cs^, F{s) EE s2 - 2Ms^ + A{sf. (107) 

Equation (fT06| now takes the form 

(108) 



dr r 

First we determine the r-coordinatc of the horizon by numerical integration, 

rn = exp ^- ^ F{s)-'^/'^A^ , (109) 

where sh = 1/(2M). (Recall that we choose ^+ , i.e. r = oo <^ s = 0, to correspond 
to r = 1.) We place the excision boundary at rmin = 0.9 rn. The ODE (jlOSp is now 
solved numerically on the interval r G [rmin, 1] with initial condition s(l) = 0. (We 
use the Python routine scipy. integrate . odeint.) 

In terms of the numerically determined function s(r), the geometry variables are 
obtained as 

fi = rs, 7r = 2Cr-3s2, N^rFisf^, X = Cs'^ -\K. (110) 
For AI = C = this reduces to Minkowski spacetime (with re [0, 1]), 

n = \K{l~r'^), 7r = 0, TV = i/i(l + r2), X = -\K. (Ill) 
For the matter fields we consider three types of initial data, 

(i) scalar field: specify initial data for (0, iji), take the Yang-Mills field to vanish, 

(ii) gravitational-sector Yang-Mills: specify initial data for {F,Df), all other fields 
vanish, 

(iii) sphaleron-sector Yang-Mills: specify initial data for {H,Dh), take F = Dp ~ 
4> ^ t/j ~ 0, solve for Dj^ and G. 

Note that in case (iii), F = Dp =0 initially but not during the evolution. In each case, 
the initial data are chosen to be a Gaussian that is approximately ingoing initially, 
e.g. for the pair (</>, "0), 

4> = Acxp(^-^-^^^^y 4>^^' + r-'4>- (112) 

For all the numerical evolutions shown here, we choose rp = 0.5 and a = 0.05. 

If we add these initial data for the matter fields to the Minkowski background 
solution (re-solving the constraints), we obtain evolutions that either disperse or 
collapse to a black hole, depending on the amplitude A. If the background geometry 
is Schwarzschild, the matter will partly accrete onto the black hole. 



6.2. Scalar field evolutions 

We begin with a scalar field evolution that disperses to flat space (figure [T]). The 
amplitude is chosen to be A = 0.6, which corresponds to an initial Bondi mass 
Mg — 0.59, already well in the non-linear regime. At late times, the field decays as a 
power law. At ^+ , </> ^ whereas at the origin and in fact at any finite radius, 
cj) ^ t~^. These results are in agreement with jTB]. The fact that the scalar field decays 
more slowly at ^+ than away from ^+ leads to the solution becoming increasingly 
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Figure 1. Scalar field dispersal. The initial Bondi mass is Mg = 0.59. At 
,J^~^ (solid line) the field decays as ~ at late times. At the origin (dashed 
line), ^ ~ 




lo' 10^ 10 



Figure 2. Scalar field accretion. The initial Bondi and apparent horizon masses 
are = 1.45, M^jj = 1, and the final masses are = M^jj = 1.44. At 
J^~^ (solid line) the field decays as </> '~ at late times. At the horizon (dashed 
fine), 4> ~ t~^. 

peaked at , similar to the formation of a boundary layer. Eventually this feature 
cannot be resolved at a fixed numerical resolution. The simulations presented here 
were run at two different resolutions in order to make sure that the plots can be trusted 
during the times shown. For the simulations shown in the plots, N = 8000 grid points 
were used. 

Next we choose initial data containing a black hole with a scalar field perturbation 
(figure 121). Initially the Bondi mass is Afg = 1.45 and the apparent horizon mass is 
M^jj = 1. Almost all of the matter falls into the black hole — the final Bondi mass, 
which agrees with the final apparent horizon mass, is = M^jj = 1.44. At late 
times we observe a power-law decay of the scalar field with the same decay exponents 
as in the evolution that dispersed to flat space, where now instead of the origin, we 
evaluate the field at the horizon. These exponents agree with those found in the test 
field approximation |24| . 

Finally we return to the case of regular initial data as in the first simulation but 
increase the amplitude such that a black hole forms in the course of the evolution 
(figure [3]). The initial Bondi mass is Mg = 2.74, and the final Bondi and apparent 
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Figure 3. Scalar field collapse. The initial Bondi mass is = 2.74 and the 
final Bondi and apparent horizon masses are Mf = A^jj = 2.71. At ,y+ (solid 
line) the field decays as </> ~ at late times. At the horizon (after it forms, 
dashed line), <f> ~ t^^. 
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Figure 4. Gravitational-sector Yang-Mills dispersal. The initial Bondi mass is 



= 0.63. At (solid lines) the fields decay as F r 
late times. At the origin (dashed lines), F ~ Dp ~ t"**. 



and Dp 



at 



horizon masses are Afg 



2.71. Again we find the same decay exponents as 



in the dispersing evolution. 

6.3. Gravitational- sector Yang-Mills evolutions 

Consider now the gravitational sector of Yang-Mills. As for the scalar field, we begin 
with an evolution that disperses to flat space (figure H]) . We find that the potential 
F decays as t~'^ at and F ^ t~'^ at the origin. These exponents agree with 
[23] . Let us also evaluate the electric field Dp. Whereas this decays at the same rate 
as F at finite radius, we observe a slower decay Dp ^ at . This may seem 
surprising at first but can be explained as follows. Consider the evolution equation 
([75]) for F, 

F = rXF' + 2XF - NDp. (113) 

Notice the radial derivative F' on the right-hand side. Suppose F' decayed at the 
same rate as F at J^"*" . Then, by continuity, F would decay at the same rate in a 
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Figure 5. Gravitational-sector Yang-Mills accretion, 
apparent horizon masses are Mi, = 1.49 and M\„ = 
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Dp 



1.45. At ^+ (solid lines) the fields decay 



The initial Bondi and 
1, and the final masses 
~ and 



at late times. At the horizon (dashed lines), F ^ Dp ^ t 




Figure 6. Gravitational-sector Yang-Mills collapse. The initial Bondi mass is 
Mg = 3.04 and the final Bondi and apparent horizon masses are M£ = M^jj = 
2.49. The final coordinate location of the apparent horizon is r\ji = 0.163. At 
,J^'^ (solid lines) F — > —1 and Dp ~ t^^ at late times. At the horizon (after it 
forms, dashed lines), F — f —37.9 ~ — and Dp ~ t~'^. 



neighbourhood of . But we know from previous work |241 123| — and our numerical 
evolutions suggest this too — that F decays faster at any point away from than 
at ^+ itself. Hence by contradiction F' must decay at a slower rate at than F. 
From (111311 we deduce that Dp must also decay at that slower rate. (This argument 
cannot determine the precise decay rate of Dp at .y^ , it merely shows that it must 
decay more slowly than F.) 

Next we consider accretion onto a black hole (figure [S]). The decay exponents 
found in this case are the same as in the dispersing evolution and agree with the test 
field approximation [24] . 

Finally we turn to collapse (figure IH]). While the electric field decays at the 
same rates as above, the behaviour of the potential F is different. At late times this 
approaches the static solution F = 2/{gr^) (recall we set g = —2). We shall see 
in section 16.51 that this is another vacuum solution of the Yang-Mills equations (in 
addition to = 0). 
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Figure 7. Sphaleron-sector Yang-Mills dispersal. The initial Bondi mass is 
= 0.75. At J'+ (solid lines) the fields decay as F ^ H ^ t~^, G ~ i"^, 
Dp ^ Dh ^ and ~ t"^. At the origin (dashed lines), F ^ H ^ G ^ 
Dp ^ Dh ^ Dl ^ 




Figure 8. Sphaleron-sector Yang-Mills accretion. The initial Bondi and apparent 
horizon masses are A/g = 1.61 and Af^^jj = 1, and the final masses are 
M| = A/|^jj = 2.49. At .y+ (solid lines), F and H approach constants, G ~ t'^ , 
Dp ^ Dh ~ and Di^ ~ At the horizon (dashed lines), F and H approach 
constants, G is frozen to its initial value G = 0, and Dp ~ Dh ~ ~ 



6.4^. Sphaleron-sector Yang-Mills evolutions 

In the last set of evolutions we eonsider initial data of type (iii) (section 16. ip that 
include the sphaleron sector of Yang-Mills. Figure [7] shows a dispersing evolution. All 
fields decay. The potentials F and H decay at the same rate as in the gravitational 
sector, i.e., F ^ H ^ t~'^ at J^+and F ^ H ^ at the origin. The "gauge 
potential" G decays as G ^ t~^ at and G ^ at the origin. The electric field 
components Dp and Dh also decay at the same rates as in the gravitational sector, 
i.e., Dp Dh ^ t-^ at J+ and Dp ^ Dh ^ t''^ at the origin. The component Dl 
has a different decay at ^+ , Dl ^ i^^, while Dl ^ t^^ at the origin as for the other 
components. 

Already in the accretion evolution (figure [5]), a different behaviour is seen for the 
connection components. Now F and H approach a nonzero static solution at late 
times. The gauge potential G is frozen to its initial value at the excision boundary, 
which for the initial data we choose is G = to within numerical roundoff. At , 
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Figure 9. Sphaleron-sector Yang-Mills collapse. The initial Bondi mass is Mg = 
3.51 and the final Bondi and apparent horizon masses are Mg = M^jj = 2.67. As 
before, solid lines refer to quantities observed at ^+ and dashed lines to quantities 
observed at the horizon after it forms. The top left panel shows that G becomes 
asymptotically time independent, G — > Go = 0.00183 at , and F and H 
perform harmonic oscillations with period T = 1717 ~ 2iT/\gGo\- In the top right 
panel, Dp and Dh decay in damped harmonic oscillations, Dj^ ~ at .y"*" and 
D]^ ~ at the horizon. The bottom panel demonstrates that all independent 
components J*", p and 5''^'''' of the energy-momentum tensor decay as t^^ at 
J^~^ and at the horizon. 



G ^ t^'^ as previously in the dispersing evolution. The decay of the eletric field is also 
the same as in the dispersing evolution. 

Finally we study sphaleron-sector collapse (figure [9]) . Recall that we freeze the 
value of G at the excision boundary once the apparent horizon forms. This results in 
the entire function G becoming time independent (but nonzero). Surprisingly however, 
the potentials F and H do not become time independent but show a sinusoidal time 
dependence with a period close to T = 2Tr/\gGo\, where Go is the asymptotic value of G 
at . This peculiar behaviour will be explained in section [6.51 A first hint at what 
is happening can be found by evaluating the components of the energy-momentum 
tensor. Figure [9] demonstrates that all of these decay at the rate at and t~^ 
at the horizon, in accordance with the electric field decaying as t~^ at J^~^ and at 
the horizon (recall that the energy-momentum tensor is quadratic in the field strength 
tensor). Hence a vacuum solution is approached. 
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6.5. The Yang-Mills vacuum 

The findings of the previous subsection suggest that we need to obtain a better 
understanding of the vacuum solutions to the Yang-Mills equations. Let us therefore 
take the field strength tensor to vanish, i.e., Dp = Dh = -Dl = Bp = Bh = Bl = 0. 
Equation ([83|) for the magnetic field implies 

rH' + 'iH = 0^ H = Ha{t)r-'^. (114) 

The linear combination x ((84|) — (|85|) yields 

rF' + 2F = {)^ F = FQity-"^. (115) 

The radial gauge condition ([89)) implies 

rG' + G = O^G = Go(<)r-i. (116) 



Since we freeze G at the inner boundary after excision, we are only interested in the 
case where Go is a constant. Substituting (|114p - (|116p into the evolution equations 
([75]) and ^ for F and H, we obtain 



Fo = -gGoHo, Ho = -Go + ^Goi^o- (117) 

The general solution to this pair of ODEs is 

Fq ^ g^^ + asm{gGot + (f)), Hn = -acos{gGot + (j)). (118) 

The constant a is not arbitrary, however. In deriving (jll7l) we used only one particular 
linear combination of (|84)) and ([85]) . Substituting the solution (|118p in either (|84l) or 
(|85p fixes a = g~^. The form of the solution (jllSp agrees well with the numerical 
results in figure [HI We see that it is Go (the value of G at J^^ ) that determines the 
period of the oscillations. If G is zero as in the final state of the accretion evolution, 
(|118p implies that Fo and Hq are constants related by 

{l-gFof + {gHo)^ = 1. (119) 

If in addition Hq = as in the gravitational-sector evolutions, this implies Fq — or 
Fo - 2/g, i.e., by 

F = 0oi F = (120) 



In the Abelian-gauge version (j90| of the gravitational-sector Yang-Mills ansatz, these 
two copies of the Yang-Mills vacuum correspond to w = ±1. 



7. Conclusions 



The purpose of this paper was to show how matter can be included in a hyperboloidal 
evolution scheme for the Einstein equations, and to implement such a scheme 
numerically in order to study power-law tails of matter fields up to future null infinity. 

We assumed that the energy-momentum tensor is tracefrce, for then the energy- 
momentum conservation equations are conformally invariant and thus regular at null 
infinity. A large class of radiative matter models are included under this assumption. 
How to deal with non-tracefree matter in the conformal setting remains an open 
question. Of course, if the matter is such that it remains bounded away from 

during the evolution then one can live with matter field equations that are singular 
at as they never need to be evaluated there. 
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We worked with a constrained ADM formulation of the Einstein equations on 
CMC shoes [3]. In that paper we showed (in the vacuum case) how the formaUy 
singular terms in the evolution equations at J^"*" can in fact be evaluated in a regular 
way. In the present paper we showed that our analysis is unaffected by the addition 
of matter (section lO)) . 

Two matter models were studied in detail, namely a conformally coupled scalar 
field and a Yang-Mills field. In both cases the energy-momentum tensor is tracefree 
and the matter equations of motion are conformally invariant. The use of conformal 
coupling instead of minimal coupling for the scalar field is essential here. We first 
derived the matter evolution equations and the source terms appearing in the Einstein 
equations without any symmetry assumptions. Subsequently we reduced the Einstein- 
matter systems to spherical symmetry in isotropic coordinates. We worked with 
the most general ansatz ([75| for the spherically symmetric Yang-Mills connection 
|331 I34[ 135] : this is more general than the purely magnetic or gravitational-sector 
ansatz that is often used. 

Our motivation to study spherically symmetric Einstein-matter evolutions on 
hyperboloidal slices arose partially from an earlier numerical study by one of the 
authors [5] that considered vacuum axisymmetric spacetimes. Whereas long-term 
stable evolutions of perturbed black holes could be achieved and the quasi-normal 
mode radiation correctly reproduced, we were unable to resolve the power-law tail of 
the gravitational field expected at late times. In the current spherically symmetric 
study a much higher numerical resolution could be used so that the tails could indeed 
be resolved. Otherwise the numerical method based on fourth-order finite differences 
is very similar to the one used in [5j . One difference is that we are now able to handle 
a regular centre as well so that gravitational collapse can be studied; in [5] we only 
considered black hole evolutions with excised interior. 

A general feature of tails in hyperboloidal evolutions is that the fields decay at 
a slower rate at ^+ than at finite radius. This means that the solution becomes 
increasingly peaked at , akin to the formation of a boundary layer. This is very 
challenging numerically — at a fixed numerical resolution the code is ultimately unable 
to resolve the solution. A possible direction for future development would be to 
introduce a time-dependent radial mapping so that an increasingly higher resolution 
can be provided near ^+ as the evolution proceeds. A similar adaptivity will be 
needed in order to study the formation of very small black holes in critical collapse. 
The current code already contains an (albeit time-independent) radial map in order 
to better resolve the steep gradients near the horizon of a black hole. 

For scalar field matter our results are consistent with evolutions of the Einstein- 
scalar field system in Bondi coordinates [18] . Our study goes further though because 
unlike Bondi coordinates, the coordinates we use can penetrate black hole horizons. 
Both for an evolution that starts out with a perturbed black hole and for an evolution 
that forms a black hole from regular initial data, we find the same decay exponents 
at J^"*" and at finite radius as in the dispersing case, which agree with results obtained 
in the test field approximation [24] . 

For evolutions of the gravitational sector of Yang-Mills, the decay rates we found 
for the Yang-Mills connection in an evolution that disperses to fiat space agree with 
evolutions in Bondi coordinates [53]. For a Schwarzschild black hole with a Yang- 
Mills perturbation we obtain the same decay rates, which in turn agree with the test 
field approximation [24] . Somewhat different behaviour is seen in an evolution that 
collapses to a black hole from regular initial data. Here the Yang-Mills connection 
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approaches a static solution that corresponds to another copy of the Yang-Mills 
vacuum. 

A fact that does not seem to have been noted before is that in all three cases 
(dispersal, accretion and collapse), the electric field has a slower decay at ^+ than 
the connection. Dp ^ as opposed to ^ t~^. At finite radius the decay rates 
agree, F ^ Dp ^ t~^. Since the electric field is the physically measurable quantity, it 
seems more natural to consider its decay rather than that of the connection. 

As far as we know we presented the first dynamical numerical evolutions that 
include the sphaleron sector of Yang-Mills, i.e., the fully general ansatz ([75)1 for 
the spherically symmetric Yang-Mills connection. In an evolution that disperses to 
Minkowski space, both potentials F and H appearing in the ansatz for the connection 
decay at the same rates as for the gravitational-sector evolutions. For an accretion 
evolution, however, both potentials approach nonzero static solutions, and in a collapse 
evolution they become sinusoidal functions of time. Nevertheless, in all cases a vacuum 
solution is approached. We were able to explain this behaviour by deriving the general 
vacuum solution to the Yang-Mills equations in our setting. 

The nontrivial structure of the Yang-Mills vacuum gives rise to interesting 
threshold behaviour [42]. The Einstein- Yang-Mills system exhibits remarkably rich 
dynamics [37[ I36j due to the existence of nontrivial asymptotically flat static solutions 
P51 133]. We hope to investigate some of these phenomena further using our 
hyperboloidal evolution code. 
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Appendix A. Conformal transformations and 3-1-1 decompositions 

In this appendix we collect a few useful identities used throughout the paper. 

Appendix A.l. Conformal transformation of the connection and curvature 
Consider a conformal transformation of the metric, 

7ab = ^^gab- (A.l) 

The following applies to any dimension n and any signature (e.g. the spacetime metric 
or the induced spatial metric on the t ~ const slices), and hence we use general indices 
a,b, . . . to indicate this. 

Let r°bc denote the Christoffel symbols of the metric connection V of 5 and V^bc 
the Christoffel symbols of the metric connection V of 7. They are related by 

r\c = - n-' (^s^n,, + s^n.b - -fbc-f'"'^dn) . (a.2) 



Hyperboloidal Einstein-matter evolution and tails for scalar and Yang-Mills fields 26 
The Ricci tensor transforms as 

Rab - Rab + {n- 2)n-^VaVb^ 

+ labf" [^-^^c^d^ -in- (A.3) 

and the Ricci scalar as 

R = ^fR + 2{n - l)VL-i''^^ b^ - n{n - l)-/°-''n^ai^^b- (A.4) 

Appendix A. 2. Physical and conformal extrinsic curvature 

The physical extrinsic curvature Kij and the unphysical extrinsic curvature Cij are 
defined as 

By translating the right-hand side of the first equation into conformal language and 
comparing with the second, we find for the traceless part 

Cl^ = nKl^^~^^~H,J,l7^'^'^'. (A.6) 

Taking instead the trace, we recover equation (|16p . 

CAn^-^{K + nc), (A.7) 

where C = 'y'''^Cij and we recall our convention K = —g'^^ Kij — const > 0. 
Appendix A.3. 3 + 1 decomposition of the scalar Hessian 

The projections of the covariant Hessian of a conformal scalar field 4> 

can be written in 3 + 1 form as 

n^n" (4)v^ (4)v^^ ^d^-N-^f^Nj^^, (A.8) 

7.^7/ ^^^V^ (4)V,0 = V,Vj^ + a,Ci,4>. (A.IO) 

Appendix A. 4-. 3 + 1 decomposition of the conformal spacetime Ricci tensor 

One of the equations of a 3 + 1 decomposition in the conformal spacetime is 

7,^7/ (^'i?^, = -Ua, ~ i-i'^'CkCji + CyC - N-^V.,VjN + R,j. (A.ll) 

Separating the trace and traceless part of Cy and using (|A.6p . we can write the 
traceless part of (jA.ll[) as 

-N-\\/'\/^Nf + 7?'"^ (A.12) 

Note the similarity of (|A.12p with formally it can be obtained from that equation 
by setting = 1, replacing K with —C and noting that the source term now refers to 
the conformal geometry. Taking instead the trace of (|A.11[) we find 

7'^" (^^i?^^ = -N-^V'V,N + R + C^ - CnC. (A.13) 

We also have the Gauss-Codazzi equation analogous to the Hamiltonian constraint 

(uni), 

^R + C^- = i? + |C2 - /i-27,fc7^,^t"J7r*'-'='. (A.14) 
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Taking linear combinations of (|A.13p and (|A.14p . wc obtain 

= -2iV-iV''V,Ar + ^ + |C2 + Ai-2^ifc7j;7r*"%*'''=' - 2£sC. (A.16) 
Finally, the Gauss-Codazzi equation analogous to the momentum constraint (|2ip 

is 

f ^n'' (^)i?^, = A^;'V,7r"-'^' + |7'-''C,,. (A.17) 
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